#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Merge & Align (Iterative)                                          #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 01/03/2007                                             #
# Last Modification: 01/03/2007                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 40
#
# MANUAL: This script allows refining the data of the selected images.
#
# MANUAL: All selected images are merged together, and subsequently aligned against the merged average. This can be done iteratively a preset number of times. Double-click onto the iteration counter after the script name, to adjust that number, if needed.
# 
#
# PUBLICATION: 3D reconstruction of two-dimensional crystals: <A HREF="http://www.ncbi.nlm.nih.gov/pubmed/26093179">Arch Biochem Biophys 581, 68-77 (2015)</A>
# PUBLICATION: 3D Reconstruction from 2D crystal image and diffraction data: <A HREF="http://dx.doi.org/10.1016/S0076-6879(10)82004-X">Methods Enzymol. 482, Chapter 4, 101-129 (2010)</A>
# PUBLICATION: 2dx - Automated 3D structure reconstruction from 2D crystal data: <A HREF="http://journals.cambridge.org/action/displayAbstract?aid=1943200">Microscopy and Microanalysis 14(Suppl. 2), 1290-1291 (2008)</A>
# PUBLICATION: 2dx_merge - Data management and merging for 2D crystal images: <A HREF="http://dx.doi.org/10.1016/j.jsb.2007.09.011">J. Struct. Biol. 160(3), 375-384 (2007)</A>
#
# DISPLAY: create_merged_dataset
# DISPLAY: MergeResolution
# DISPLAY: zstarwin
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: merge_res_limit
# DISPLAY: tempkeep
# DISPLAY: realcell
# DISPLAY: realang
# DISPLAY: ALAT
# DISPLAY: MergeStepSize
# DISPLAY: IBOXPHS
# DISPLAY: SYM
# DISPLAY: avrgamphsNUMBER
# DISPLAY: avrgamphsRESOL
# DISPLAY: MergeIQMAX
# DISPLAY: Merge_Reference_IQMAX
# DISPLAY: MergeHKMAX
# DISPLAY: Merge_Reference_HKMAX
# DISPLAY: AMP_Scale_Factor
# DISPLAY: ILIST
# DISPLAY: refbeamtilt
# DISPLAY: reftiltgeo
# DISPLAY: ITAXASTEP
# DISPLAY: RTAXASIZE
# DISPLAY: ITANGLSTEP
# DISPLAY: RTANGLSIZE
# DISPLAY: MergeAK
# DISPLAY: MergeIWF
# DISPLAY: MergeIWP
# DISPLAY: merge_ref_num
# DISPLAY: merge_comment_1
# DISPLAY: merge_comment_2
# DISPLAY: merge_comment_3
# DISPLAY: merge_comment_4
# DISPLAY: merge_comment_5
# DISPLAY: merge_comment_6
# DISPLAY: merge_comment_7
# DISPLAY: merge_comment_8
# DISPLAY: merge_comment_9
# DISPLAY: merge_refine_iterations
# DISPLAY: max_amp_correction
# DISPLAY: make_reference
# DISPLAY: merge_data_type
# DISPLAY: plotres_rings
# DISPLAY: resolutionplot_RESMAX
# DISPLAY: resolutionplot_bins
# DISPLAY: RFACAMP
# DISPLAY: Thread_Number
# DISPLAY: RESMIN
# DISPLAY: RESMAX
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
#
set ALAT = ""
set AMP_Scale_Factor = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set Calc_from_zstarrange = ""
set create_merged_dataset = ""
set IBOXPHS = ""
set ILIST = ""
set ITANGLSTEP = ""
set ITAXASTEP = ""
set lattice = ""
set max_amp_correction = ""
set merge_data_type = ""
set merge_ref_num = ""
set merge_reference = ""
set Merge_Reference_HKMAX = ""
set Merge_Reference_IQMAX = ""
set merge_refine_iterations = ""
set merge_res_limit = ""
set MergeAK = ""
set MergeHKMAX = ""
set MergeIQMAX = ""
set MergeIWF = ""
set MergeIWP = ""
set MergeResolution = ""
set MergeStepSize = ""
set num_reflections_fitted = ""
set num_reflections_FOM50 = ""
set plotres_rings = ""
set realang = ""
set realcell = ""
set refbeamtilt = ""
set reftiltgeo = ""
set RESMAX = ""
set RESMIN = ""
set resolutionplot_bins = ""
set resolutionplot_RESMAX = ""
set RFACAMP = ""
set RTANGLSIZE = ""
set RTAXASIZE = ""
set SYM = ""
set tempkeep = ""
set Thread_Number = ""
set zstarwin = ""
set ctfrev = ""
#
#$end_vars
#
setenv OMP_NUM_THREADS ${Thread_Number}
#
#
set scriptname = 2dx_mergeRefine2D
\rm -f LOGS/${scriptname}.results
#
set merge_modus="2D"
set merge_modus = "2D"
set refbeamtilt = n
set reftiltgeo = n
set merge_reference = 0
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set dirfile = "2dx_merge_dirfile.dat"
set dirnum = `cat ${dirfile} | wc -l`
set maxthread = `echo ${Thread_Number} ${dirnum} | awk '{if ($1<$2/2) { s = $1 } else { s = int($2 / 2) }} END { print s }'`
if ( ${maxthread} == "0" ) then
  set maxthread = 1
endif
#
set number = 1
if ( ${ILIST} == "n" ) then
  set IVERBOSE = 1
else
  set IVERBOSE = 6
endif
#
#
# The following is to make sure that for the next "Import Images", the default is correctly initialized.
set initialization_reset = "y"
set initialization_executable = "y"
echo "set initialization_reset = ${initialization_reset}" >> LOGS/${scriptname}.results
echo "set initialization_executable = ${initialization_executable}" >> LOGS/${scriptname}.results
#
if ( ${merge_ref_num} != "0" ) then
  ${proc_2dx}/linblock "ERROR: Iterative refinement does not make sense with"
  ${proc_2dx}/linblock "a saved static reference."
  ${proc_2dx}/linblock "Use last merge result as reference (register 0),"
  ${proc_2dx}/linblock "or use the REFINE script instead."
  ${proc_2dx}/protest "Aborting."
endif
#
echo MergeAK = ${MergeAK}
echo MergeIWF = ${MergeIWF}
echo MergeIWP = ${MergeIWP}
#
echo "<<@progress: 1>>"
#
#################################################################################
${proc_2dx}/linblock "Verifying some parameters"
#################################################################################
#
if ( ${MergeHKMAX}x == 'x' ) then
  set MergeHKMAX = '20'
  ${proc_2dx}/linblock "ERROR: correcting MergeHKMAX to ${MergeHKMAX}"
  echo "set MergeHKMAX = ${MergeHKMAX}" >> LOGS/${scriptname}.results
endif
#
if ( `echo ${RESMIN} ${RESMAX} | awk '{ if ( $1 < $2 ) { s = 1 } else { s = 0 }} END { print s }'` == 1 ) then
  set oldval = ${RESMIN}
  set RESMIN = ${RESMAX}
  set RESMAX = ${oldval}
  ${proc_2dx}/linblock "ERROR: exchanging RESMIN and RESMAX, to RESMIN=${RESMIN}, and RESMAX=${RESMAX}"
  echo "set RESMIN = ${RESMIN}" >> LOGS/${scriptname}.results
  echo "set RESMAX = ${RESMAX}" >> LOGS/${scriptname}.results
endif
#
if ( ${AMP_Scale_Factor}x == 'x' ) then
  set AMP_Scale_Factor = 0.001
  echo "set AMP_Scale_Factor = ${AMP_Scale_Factor}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: AMP_Scale_Factor corrected to ${AMP_Scale_Factor}"
endif
# 
if ( ${MergeAK}x == 'x' ) then
  set MergeAK = "10"
  echo "set MergeAK = ${MergeAK}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: MergeAK corrected to ${MergeAK}"
endif
#
if ( ${MergeIWF}x == 'x' ) then
  set MergeIWF = "0"
  echo "set MergeIWF = ${MergeIWF}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: MergeIWF corrected to ${MergeIWF}"
endif
#
if ( ${MergeIWP}x == 'x' ) then
  set MergeIWP = "1"
  echo "set MergeIWP = ${MergeIWP}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: MergeIWP corrected to ${MergeIWP}"
endif
#
if ( ${ILIST}x == 'x' ) then
  set ILIST = "n"
  echo "set ILIST = ${ILIST}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: ILIST corrected to ${ILIST}"
endif
#
if ( ${ILIST} == "n" ) then
  set ILIST_VAL = 0
else
  set ILIST_VAL = 1
endif
#
if(${MergeIWF} == "0")then
  set MergeIWF_VAL = -1
endif
if(${MergeIWF} == "1")then
  set MergeIWF_VAL = 0
endif
if(${MergeIWF} == "2")then
  set MergeIWF_VAL = 1
endif
if(${MergeIWP} == "0")then
  set MergeIWP_VAL = -2
endif
if(${MergeIWP} == "1")then
  set MergeIWP_VAL = -1
endif
if(${MergeIWP} == "2")then
  set MergeIWP_VAL = 0
endif
if(${MergeIWP} == "3")then
  set MergeIWP_VAL = 1
endif
echo MergeAK = ${MergeAK}
echo MergeIWF_VAL = ${MergeIWF_VAL}
echo MergeIWP_VAL = ${MergeIWP_VAL}
#
# This translates the list of directories to work on into one single long line:
cat 2dx_merge_dirfile.dat | tr "\n" " " > SCRATCH/2dx_merge_dirfile_oneline.dat
set dirlist = "`cat SCRATCH/2dx_merge_dirfile_oneline.dat`"
#
# This memorizes the current merge directory under the variable "olddir":
set olddir = $PWD
#
echo "The current working directory is" ${olddir}
#
echo "<<@progress: 5>>"
#
#############################################################################
${proc_2dx}/linblock "Sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
############################################################################# 
${proc_2dx}/lin "2dx_merge_makedirs - to create all required subdirectories"
#############################################################################
#
source ${proc_2dx}/2dx_merge_makedirs
#
echo "<<@progress: 10>"
#
if ( ! -e APH/merge.aph ) then
    set create_merged_dataset = "y"
endif
#
set NPRG = 1
if ( ${merge_reference} == '0' ) then
  # Use merge.aph files
  if ( ! -e APH/merge.aph ) then
    ${proc_2dx}/linblock "ERROR: APH/merge.aph not found. No reference for refinement available."
    ${proc_2dx}/protest "Run Merging script first."
  endif
  set reference_file="APH/merge.aph"
endif 
if ( ${merge_reference} == '1' ) then
  # Use interpolated lattice lines merge2Dref_MRClefthanded.mtz
  if ( ! -e merge2Dref_MRClefthanded.mtz ) then
    ${proc_2dx}/linblock "ERROR: merge2Dref_MRClefthanded.mtz not found. No reference for refinement available."
    ${proc_2dx}/protest "Run Merging script first."
  endif
  set reference_file="merge2Dref_MRClefthanded.mtz"
  set NPRG = 3
endif
#
set itnum = 0
#
while ( ${itnum} < ${merge_refine_iterations} ) 
  #
  set itnum_p1 = ${itnum}
  @ itnum_p1 += 1
  echo ":: "
  echo "::         ITERATION ${itnum_p1} of ${merge_refine_iterations}"
  echo ":: "
  #
  if ( ${merge_ref_num} == "0" && ${create_merged_dataset} == "y" ) then
    #
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_merge.com"
    #############################################################################  
    #
    sleep 2
    source ${proc_2dx}/2dx_origtilt_merge.com
    #
  endif
  #
  #############################################################################
  ${proc_2dx}/linblock "Sourcing 2dx_origtilt_refine.com"
  #############################################################################  
  #
  echo "PWD = "$PWD
  source ${proc_2dx}/2dx_origtilt_refine.com
  #
  echo "<<@evaluate>>"
  #
  @ itnum += 1
end
#
${proc_2dx}/linblock "Now, press the REFRESH DISPLAY VIEW button in top menu to update the project library."
#############################################################################
${proc_2dx}/linblock "2dx_refine normal end."
#############################################################################
#
echo "<<@progress: 100>>"
#
# 
exit
#
